clc,clear;
close all;
%导入相关数据
data=xlsread("random_coordinates.xlsx");
% 城市坐标
cities=data*100;
% 城市数量
nCities = length(cities);
%计算距离矩阵
dist = zeros(nCities, nCities);
for i = 1:nCities
    for j = 1:nCities
        dist(i, j) = sqrt(sum((cities(i, :) - cities(j, :)).^2));
    end
end
maxiter=1000; %最大迭代次数
startcity=5;    %初始城市
temperature=1000;    %初始温度
ftemperature=0.01;    %最终温度
rate=0.99;         %温度变化率
energy_d=8;        %能量阈值
[bestx,bestf,bd]=improved_sa_tsp(nCities,dist,startcity,temperature,rate,ftemperature, maxiter, energy_d);
disp("Best Solvution");
disp(bestx);
disp("Best Function");
disp(bestf);
tusi(bd,nCities,cities,bestx);

function [bestr, bestf,bd] = improved_sa_tsp(ncities, dist, startcity, temperature, rate, ftemp, max_iter, energy_d)
    % 初始化路径，将指定的起始城市固定在路径的首位
    remaining_cities = setdiff(1:ncities, startcity);
    croute = [startcity, remaining_cities(randperm(length(remaining_cities)))];
    cdistance = caldistance(croute, dist);
    
    % 设置初始最优路径
    bestr = croute;
    bestf = cdistance;
    
    % 初始温度和能量
    temp = temperature;
    energy = 0;
    j=1;
    while temp > ftemp
        for i = 1:max_iter
            % 生成新路径（通过交换、移位、倒置），但固定首个城市
            new_route = perr(croute);
            new_distance = caldistance(new_route, dist);
            
            % 计算增量
            delta = new_distance - cdistance;
            
            % 根据SA准则决定是否接受新路径
            if delta < 0 || exp(-delta / temp) > rand
                croute = new_route;
                cdistance = new_distance;
            end
            
            % 更新最优路径
            if cdistance < bestf
                bestr = croute;
                bestf = cdistance;
            end
            
            % 更新能量
            if delta < 0
                energy = energy + 1;
            end
        end
        
        % 跳跃式降温
        if energy > energy_d
            temp = temp * (rate ^ 2);  % 快速降温
        else
            temp = temp * rate;  % 正常降温
        end
        
        % 重置能量
        energy = 0;
        bd(j)=bestf;
        j=j+1;
    end
    
    % 二重退火：重启退火过程
    temp = temperature;
    croute = bestr;
    cdistance = bestf;
    energy = 0;
    while temp > ftemp
        for i = 1:max_iter
            % 生成新路径，但固定首个城市
            new_route = perr(croute);
            new_distance = caldistance(new_route, dist);
            
            % 计算增量
            delta = new_distance - cdistance;
            
            % 根据SA准则决定是否接受新路径
            if delta < 0 || exp(-delta / temp) > rand
                croute = new_route;
                cdistance = new_distance;
            end
            
            % 更新最优路径
            if cdistance < bestf
                bestr = croute;
                bestf = cdistance;
            end
            
            % 更新能量
            if delta < 0
                energy = energy + 1;
            end
        end
        
        % 跳跃式降温
        if energy > energy_d
            temp = temp * (rate ^ 2);  % 快速降温
        else
            temp = temp * rate;  % 正常降温
        end
        
        % 重置能量
        energy = 0;
        bd(j)=bestf;
        j=j+1;
    end
    %完整路径
    bestr=[bestr,startcity];
end

function tdist = caldistance(route, dist)
    % 计算给定路径的总距离，并确保路径返回起点
    tdist = 0;
    for i = 1:length(route)-1
        tdist = tdist + dist(route(i), route(i+1));
    end
    tdist = tdist + dist(route(end), route(1));  % 封闭路径回到起点
end

function newroute = perr(route)
    % 路径变异：交换、移位、倒置，但固定首个城市
    method = randi(3); % 随机选择变异方式
    start_city = route(1);
    recities = route(2:end);
    
    switch method
        case 1 % 交换
            idx = randperm(length(recities), 2);
            new_recities = recities;
            new_recities(idx(1)) = recities(idx(2));
            new_recities(idx(2)) = recities(idx(1));
        case 2 % 移位
            idx = sort(randperm(length(recities), 2));
            new_recities = recities;
            segment = new_recities(idx(1):idx(2));
            new_recities = [new_recities(1:idx(1)-1), segment(end), segment(1:end-1), new_recities(idx(2)+1:end)];
        case 3 % 倒置
            idx = sort(randperm(length(recities), 2));
            new_recities = recities;
            new_recities(idx(1):idx(2)) = recities(idx(2):-1:idx(1));
    end
    
    newroute = [start_city, new_recities];
end

%% 作图

function ktu=tusi(l,nCities,cities,bP)
l=[l(2:end)];
n=length(l);
count=1:n;
x=zeros(nCities+1,1);
y=zeros(nCities+1,1);
for i=1:nCities
    x(i)=[cities(bP(i),1)];
    y(i)=[cities(bP(i),2)];
end
x(i+1)=cities(bP(1),1);
y(i+1)=cities(bP(1),2);

figure(1)
plot(count,l);
xlabel('迭代次数');
ylabel('最优值');
title('模拟迭代图');
hold on
figure(2)
plot(x,y,'-ro', 'LineWidth', 2, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
xlabel('X 坐标');
ylabel('Y 坐标');
title('路径图');
grid on; % 打开网格
end

